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摘要 当 史 瓦 西 黑洞 周围 存在 渐 近 均匀 的 外 部 磁场 时 , 描述 带电 粒子 在 史 瓦 西 黑洞 附近 运动 的 哈密 顿 系统 会 变 
为 不 可 积 系统 . 类 似 于 这 样 的 相对 论 哈密 顿 系 统 不 存在 有 显 式 分 析 解 的 2 部 分 分 离 形式 , 给 显 式 辛 算法 的 构建 和 
点 用 带 来 困难 . 近 一 年 以 来 的 系列 工作 提出 将 相对 论 哈密 顿 系统 分 解 为 具有 显 式 分 析 解 的 2 个 以 上 分 离 部 分 形 
式 , 成 功 解决 了 许多 相对 论 时 空 构建 显 式 辛 算法 的 难题 . 最 近 的 工作 回答 了 哈密 顿 系 统 显 式 可 积分 离 数 目 对 长 
期 数值 积分 精度 有 何 影 响 、 哪 种 显 式 辛 算法 有 最 佳 长 期 数值 性 能 这 两 个 问题 , 指出 哈密 顿 有 最 小 可 积分 离 数目 
即 3 部 分 分 裂解 形式 并 且 应 用 于 优化 的 4 阶 分 段 龙 格 库 塔 显 式 辛 算法 可 取得 最 好 精度 . 由 此 选择 上 述 数值 积分 方 


i 


法 并 利用 庞 加 莱 截 面 、 最 大 李 雅 普 诺 夫 指数 和 快速 李 雅 普 诺 夫 指标 研究 在 磁化 史 瓦 西 黑洞 附近 运动 的 带电 粒子 
轨道 动力 学 . 结果 显示 : 针对 某 特 定 的 粒子 能 量 和 角 动 量 , 较 小 的 外 部 磁场 很 难 形成 混沌 轨道 ; 较 大 的 正 磁场 参 
数 容易 使 轨道 产生 混沌 , 并 且 随 着 磁场 的 增 大 , 轨道 的 混沌 程度 也 随 之 加 强 ; 粒子 能 量 适当 变 大 也 可 以 加 剧 混 沌 
程度 , 但 负 磁 场 参 数 和 粒子 角 动 量变 大 都 会 减弱 混沌 . 
关键 词 天 体力 学 , 黑洞 , 磁场 , 混沌 , 算法 
中 图 分 类 号 : P138; ”文献 标识 码 : A 

1 引言 到 的 黑洞 都 有 等 离子 体 吸 积 盘 , 可 以 在 黑洞 外 产生 


近 几 年 来 ,引力 波 上 3 的 数 十 起 成 功 探测 以 及 磁场 . 这 样 的 磁场 一 般 比 较 弱 小 , 虽然 对 黑洞 时 空 
M87 星 系 中 心 超大 质量 黑洞 (M87*) 照 片 8 扫 和 银河 没有 什么 影响 , 但 对 蓓 质 比 (电荷 与 质量 之 比 ) 较 大 
系 中 心 黑洞 (SgrA*) 照 片 5-9 陆 续 证 实 了 爱 因 斯 坦 的 粒子 产生 洛 伦 兹 力作 用 对 粒子 运动 影响 不 可 忽 
So 广义 相对 论 引力 理论 关于 黑洞 存在 的 预言 . 黑洞 物 ME. 当 史 瓦 西 时 空 有 外 部 电磁 场 浸入 时 , 哈密 
> 里 性 质 和 引力 检验 已 成 为 研究 热点 问题 . 顿 雅 可 比方 程 很 可 能 会 变 为 不 可 分 离 变 量 的 情形 ， 
史 瓦 西 黑洞 是 爱 因 斯 坦 场 方程 的 一 个 解 . 由 于 即 不 存在 类 似 的 卡特 常数 , 从 而 相应 的 哈密 顿 系统 
史 瓦 西 黑洞 存在 4 个 独立 的 运动 常数 , 即 粒 子 能 量 、 变 为 不 可 积 系统 . 带电 粒子 在 磁化 黑洞 周围 的 运动 


角 动 量 、 粒 子 的 4 速度 关系 和 类 似 的 卡特 常数 ( 亦 称 定 条 件 下 会 产生 混沌 现象 中 ?9. 判定 试验 粒 
KAMAE, 所 以 是 可 积 系统 . 尽管 该 系统 理 子 轨道 的 混沌 性 方法 有 庞 加 莱 截 面 、 最 大 李 雅 普 
论 上 存在 解析 解 , 但 分 析 解 很 难 用 初等 函数 表示 ， 诺 夫 指数 2 和 快速 李 雅 普 诺 夫 指标 多 等 . 

仅仅 可 以 用 椭圆 积分 形式 表达 而 已 . 大 多 数 可 观测 数值 方法 是 研究 混沌 系统 最 常用 的 方法 . 对 于 


2022-04-08 收 到 原稿 , 2022-05-31 收 到 修改 稿 
* 国 家 自然 科学 基金 项 目 (11973020) 资 助 
1M130120101@sues.edu.cn 


Chinax iv ERAT! 


64 卷 天 文 学 报 4 期 


哈密 顿 系统 长 期 数值 积分 方法 来 说 ,保持 哈密 顿 PSS 

相 流 的 辛 积分 器 3 自然 成 为 优先 选择 ,可 使 能 量 在 球 坐标 (t,7,0, wp) 中 ,表示 坐标 时 间 , +r、0、 
A REECE. 辛 算法 分 为 显 式 辛 算 。 ”表示 3 维 空间 , 类 似 于 空间 直角 坐标 系 中 的 zx、 人 外 
法 B3-24 和 隐 式 辛 算法 9; 显 隐 混 合算 法 R56 纪 也 常 ”> 描述 一 个 带电 荷 g 的 试验 粒子 在 磁化 史 瓦 西 黑洞 
被 使 用 . 显 式 算法 比 隐 式 方法 具有 计算 效率 优势 ， ”附近 运动 的 动力 学 哈密 顿 函数 方程 为 : 
相对 论 时 空 对 应 的 哈密 顿 一 般 没 有 显 式 分 析 解 的 Eag 

两 部 分 分 离 形式 , 导致 显 式 辛 算法 的 构建 和 应 用 有 Am mA 
极 大 困难 . 求解 这 些 相对 论 时 空间 题 , 如 采用 辛 算 ”其 中 , g 是 粒子 的 电荷 , PF 是 由 一 组 标准 哈密 顿 正则 
法 的 话 , 过 去 只 能 是 隐 式 [5 或 者 显 隐 混 合 2631 的 。 ”方程 i = 级 确定 的 广义 动量 的 . p, p RR 
ER. 当然 , 结合 中 点 置换 的 扩大 相 空 间 显 式 类 辛 BANJE, 4A,,、A, 表 示 电 磁 四 向 量 势 , 上 下 标 从 、 
方法 04 32-371] 可 以 适合 用 于 求解 不 能 变量 分 离 的 相 wv 表示 不 同 坐 标 符号 , 代表 球 坐 标 (t,7,0, yp) 中 任意 
对 论 时 空间 题 . 近 一 年 以 来 伍 均 教授 团队 在 天 体 物 。 ”化 标 符号 , 史 瓦 西 黑 洞 度 规 的 非 零 逆 变 分 量 9 为 : 
音 杂 志 等 期 刊 上 发 表 了 系列 工作 07, 20, 38-44) pi 2\-! 

将 相对 论 时 空 对 应 的 哈密 顿 或 时 间 变 换 的 哈密 屯 =1ou=-(1-=) ， 
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Z 分 解 为 具有 显 式 分 析 解 的 ?个 以 上 分 离 部 分 形式 ， Pop 
O ”成功 解 决 了 许多 相对 论 时 空 构建 显 式 辛 算法 的 难 r 
D 题 这 一 思路 还 可 以 应 用 于 Yoshida 高 阶 显 式 六 算 g = 100 =, 9°? = 1/90 = =a. (2) 
SMO 问题 是 哈密 顿 系统 显 式 可 积分 离 数 目 是 。 本 
O 否 影响 长 期 数值 积分 精度 ， 还 有 一 个 问题 是 何 种 这 里 gu gers go> doo FE SALTER IEE SEN ME 
O 显 式 辛 算法 有 最 佳 长 期 数值 性 能 .作者 最 近 发 表 ”分 最 a ee 
Rik LA NINET EAA Aid RE EAL AR IN TAY SRL. 两 个 常数 
NG, SEMA RTBU AE Rama PAINAA: 
© A, JA nin WA AX ATS Hb ID TS 4S 
A 解 形式 并 且 应 用 于 优化 的 4 阶 分 段 龙 格 库 塔 显 式 六 en = (a 2 ) P, (3) 
> 7% (Partitioned-Runge-Kutta, PRK,4)!47) ay Rx i 
最 好 精度 . 注意 这 种 4 阶 PRK 方 法 与 同 阶 Yoshida 算 Pe = Sep tay =r sin? 0 (9+ 22) = 工 
法 相 比 , 包含 更 多 领 外 的 时 间 系数 和 子 哈密 顿 解 的 4) 
组 合 ， 故 优化 后 的 PRK 算 法 在 计算 上 比 Yoshida 算 a 、_ _ 
法 更 耗 时 . 但 额外 增加 的 时 间 不 太 多 , 并 且 拥 有 更 Po Pp AATEC AM LIN Lalit, BEAN 
小 的 截断 误差, 因而 这 种 方法 值得 推荐 应 用 . 磁场 强度 ，E 是 粒 了 的 恒定 能 量 , 工 是 粒子 的 恒定 
本 文 主要 目的 是 沿 着 近期 工作 [9 的 思路 , RATE. AP RRS IE SOR RR 
]4 阶 优化 的 PRK 显 式 辛 方法 应 用 于 具有 显 式 分 析 MO ER: 
解 的 3 部 分 哈密 顿 分 解 形 式 并 结合 庞 加 菜 截面 、 最 Ay = Č gpp = Žr? sin? 8. (5) 
大 李 雅 普 诺 夫 指数 和 快速 李 雅 普 诺 夫 指标 研究 在 


人 将 (3) (5) 式 代入 () 式 ,磁化 中 瓦 西 哈密 顿 可 以 简 
磁化 史 瓦 西 黑洞 附近 运动 的 带电 粒子 轨道 动力 学 
i A 


st p 
2 ”物理 模型 CEEE 
本 节 先 给 出 描述 带电 粒子 在 附带 渐 近 均匀 磁 2 r 2 i 
场 的 史 瓦 西 黑洞 周围 运动 的 动力 学 模型 , 再 给 出 适 De | 1 (: wei) = 

合 于 求解 磁化 史 瓦 西 时 空 哈密 顿 系统 的 辛 算法 . 2r? ' 2r?sin? 0 2 
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He 
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Dr、26 分 别 表示 在 "与 方向 上 的 广义 动量 ， 这 


里 6 = qB. (6) 式 实际 是 经 过 无 量 纲 化 处 理 后 所 得 
到 的 . 光速 c 和 引力 常数 G 单 位 化 处 理 : c= G = 1, 


量 纲 处 理 就 是 以 黑洞 质量 M 和 粒子 质量 mm 来 进 
行 标 度 变 换 : t > tM,r > rM,B > B/M,E > 
mE, p, > mp,,L > mML,p > mMpg,q 一 
mg, H 一 72 万 .其 中 mm 为 试验 粒子 质量 ，M 引 力 
源 质量 变 为 几何 化 单位 , 即 : M = 1. 

哈密 顿 函 数 (6) 式 除 含 两 个 常数 媚 和 元 外 ,还 有 
哈密 顿 本 身 总 是 一 个 守恒 量 . 对 于 类 时 时 空 , 这 个 
常数 是 : 


E. 
=) 


1 


H=-3. (7) 


当 史 瓦 西 黑洞 外 无 渐 近 均匀 磁场 时 , 哈密 顿 函 
数 (6) 式 的 第 4 个 运动 积分 存在 , 故而 是 严格 可 积 的 ; 
但 当 存在 渐 近 均匀 磁场 时 ,第 4 个 运动 积分 不 存在 ， 
从 而 导致 系统 不 可 积 , 进而 可 能 会 产生 混沌 . 


2.2 ”算法 实现 

文献 [39] 将 哈密 顿 函 数 (6) 式 分 解 为 4 个 显 式 可 
积 求解 部 分 构建 显 式 辛 算法 . 作者 近期 的 工作 [的 指 
出 哈密 顿 函 数 (6) 式 还 可 分 解 为 3、5 个 显 式 可 积 求 
解 部 分 来 构建 显 式 辛 算 法 ,发 现 3 部 分 分 解 方法 精 


度 最 高 . 按照 这 一 思路 , 3 个 部 分 裂解 方法 如 下 : 
H = Hı + Hə + H3, 


(8) 


ar 
1 ( qBr? sin? *) ? (9) 
2r? sin? 0 2 , 
H, = ip + La (10) 
Hy =—"p2. (11) 


很 明显 , 3S ERPS tee LA. Has Ath EAR 
可 以 写 出 解析 解 ， 其 解析 解 是 时 间 t 的 显 式 函数 . 
SRA ZEB Oy EEF RYE MAH. H, H. 


设 h 为 积分 步 长 , MMS R H GE PRA — We A 
Ff: 
SIH = Hi x HIP x HE, (12) 


电 测 试 粒子 在 磁化 史 瓦 西 黑洞 中 的 混沌 运动 
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S1H = HF x HE x HI. (13) 


其 中 S1H* 是 S1H 的 共 轿 辛 算 子 . 于 是 , 这 2 个 一 阶 辛 
算 子 的 积 是 一 个 对 称 组 合 , 也 就 是 二 阶 显 式 辛 方法 : 


3 个 2 阶 对 称 乘积 可 以 构造 4 阶 Yoshida 显 式 辛 
积分 器 上 5]: 


h 


S2H = S1H (4) x S1H* ( (14) 


S4H = S2H(c,h) x S2H(cgh) x S2H(cih), (15) 
AP: cl =1/(2 — 243), o =1 — 2c. 进一步 可 


以 建立 优化 的 4 阶 分 段 龙 格 库 塔 PRK64 显 式 辛 算 
法 [4 47], 


PRK64 = S1Ha,, x SIH}, x S1Ha,, x SIHG, x 
S1Has x SIH, x S1Has x SIH, x 
S1Ha, x SIH. x S1Ha, x SLE... 


10 


(16) 
从 第 (n 一 1) 步 到 第 n 步 , 具体 离散 差分 方案 如 
下 : 
S1H 
HE :rl = [(r2 4 site) eal 
PL? = Prm—il(Paa hppa ral, 
Hi? :el = C08 On_1/7r™* + pl sin On_1/D9n—1 5 
ez = sinO,_1/r? — pl cos On_1/Pen—1; 
fı = arctan 2(e2, €1) , fo = tan(@,_1 — fı), 
072 = fi + arctan[(e? + e3)hpen—1 + fol, 
Dr? = Po.n—1(e€1 SÌN Ôn—1 — €2 COS Ôn—1), 
rẹ? = 1/(e1 cos0 1 十 ezsin0 1). 
Hi : pit = p + ALLL — B(r#?)? sin? 6%? /2]? 
oy sin? 7? |= 
{B[L — B(r™)? sin? 97 /9)} /r 2 
EB? /(r? — 2)°},, 
pe = Pon—1 + h[—(cos07*{L — [B(r#?)? 
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sin? 6 /2]}) /[(r#2)? sin? 2] 
{8 cos 0™®[L — B(r™#)? 
sin? 0 /2]}/ sin 0™?]. 
S1H* 

HE spt = pen thE Bras 
sin? 6,_1/2)*/(r3_, sin? 0,,_1)— 
[B(L — Bra- sin? On—1/2)]/tr1— 
B /(rn-1 — 2)"), 

Po" = Pon—1 + h{—[cos O,_1(L — (Brn 

sin? 0,-1/2)")]/(r3_1 sin? 0,_1)— 
[8 cos On—1(L — Bra_i 
sin? 0,-1/2)]/ sin ,_1] . 


*Hy * Ay 


qi, ef = COS On_1/Tn-1 + pe? sin On_1/p5°°* 5 
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cont) Jay ， 
fi), 
gm 一 ft + arctan|(e%? + e% ape” + fol, 


e3 = sin On-1/Tn-1 — Py 
fi = arctan 2(e3, eï) , f3 = tan(@n-1 — 
f= on (ej sin 8,1 — es cosO,_1), 
r*H2 一 1/(e} cos0 1 十 ezsinb0 1). 

pris — {[(rt#2)? = Shpt te]? [pt Hay 1/3 ， 
pits pie (er)? aps] (PPS, 


*H3 è 
Hrs. 


式 中 计算 出 的 ai, Qo,…… , ai2 列 举 在 表 1 中 . 所 谓 优 
化 算法 意味 着 系数 a;、b; 之 间 的 自由 系数 在 5 阶 ( 优 
化 的 4 阶 显 式 辛 算法 ) 截 断 误差 项 中 系数 平方 和 最 
小 化 . 与 非 优 化 算法 相 比 ， 优 化 算法 可 以 减 小 
截断 误差 . (12)-(16) 式 都 是 求解 哈密 顿 函 数 (6) 式 
或 (8) 式 的 显 式 辛 算法 . 


#1 PRK64 算 法 相关 系数 
Table 1 Correlation coefficient of PRK64 integrators 


Integrators 


PRKe4 


Q1 = Q12 = 0.079203696431196 


Coefficient 


Qs5 = Ag = 0.324648188689706 


a3 = Q10 = 0.222861495867608 


a2 = Q11 = 0.130311410182166 
a4 = ag = —0.366713268047426 


Qe = a7 = 0.109688477876750 


3 数值 分 析 


辛 算法 PRK64. 可 以 看 到 两 条 曲线 都 不 随时 间 的 
增 大 而 增 大 ， 在 长 期 积分 中 保持 平稳 ， 显 式 辛 


先 评估 算法 优 劣 , 再 挑选 好 的 算法 并 结合 混 六 
指标 来 研究 带电 粒子 轨道 的 动力 学 , 特别 探讨 3 个 
动力 学 参数 对 轨道 混沌 的 影响 . 

3.1 ”算法 对 比 

取 时 间 步 长 为 h = 1 以 及 参数 为 磁场 参数 8 = 
8.9x1074ł, MT fE ÆE = 0.995 和 和 角 动量 L = 4.6. 
初始 条 件 为 p = 0、0 = 2/2. 一 旦 初始 半径 7 给 
定 ， po 初 值 取 正 值 并 且 由 (7) 式 确定 . 图 1 为 两 
种 4 阶 算法 的 哈密 顿 误差 图 ， 图 1 (a) 为 两 种 4 阶 
算法 积分 初始 半径 r = 11 的 轨道 所 得 哈密 顿 误 
#ŽAH = H + 1/28, 其 中 红色 的 曲线 表示 Yoshida 


[eS 


算法 $4 的 哈密 顿 误差 阶 数 值 最 后 稳定 在 10 的 -9 阶 
到 -8 阶 次 之 间 ， 而 优化 后 的 显 式 辛 算法 PRK64 算 
法 的 哈密 顿 误 差 阶 数值 稳定 在 10 的 -13 阶 到 -12 阶 
次 之 间 ， 比 S4 算 法 精度 高 4 个 数量 级 . 图 1 (b) 为 初 
始 半径 > = 88 的 误差 图 . 两 种 4 阶 算法 的 精度 非常 
高 . 大 轨道 初始 半径 比 小 轨道 初始 半径 精度 高 的 
原因 在 于 前 者 平均 轨道 周期 要 小 于 后 者 ， 但 积分 
到 7 = 104 后 两 条 曲线 全 部 上 移 , 原因 是 截断 误差 
远 小 于 计算 机 的 舍 入 误差 导致 的 . 即使 误差 有 长 期 
增长 趋势 ,但 精度 仍 可 达到 10 的 -13 阶 次 左右 . 还 
要 注意 到 优化 后 的 显 式 辛 算法 需要 使 用 更 多 的 计 


4 阶 显 式 辛 算法 S4， 绿 色 曲 线 表示 PRK 4 阶 显 式 


算 时 间 . 表 2 中 列 出 两 种 算法 的 中 央 处 理 器 (Central 
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Processing Unit, CPU) 时 间 , 表明 优化 后 的 PRK 方 
法 计算 时 间 比 S4 算 法 稍 长 一 些 , 但 总 的 来 说 , 计算 
时 间 依 旧 很 短 , 最 长 CPU 时 间 不 超过 3 min. 因此 ， 
考虑 到 精度 和 效率 ，PRK64 算 法 值得 在 实 算 中 使 


(b)r= 88 


Ig lAH| 


41 两 种 4 阶 算法 的 哈密 顿 误差 图 . 磁场 参数 8 = 8.9x1074, 粒子 能 
HE = 0.995, 角 动 量 L = 4.6. 


Fig.1 Hamiltonian errors for two fourth-order integrators. 
Magnetic parameter 6 = 8.9x10~“, particle energy 


E = 0.995, angular momentum L = 4.6. 


表 2 图 1 中 两 种 算法 的 CPU 时 间 
Table 2 CPU times for the two integrators in 


Fig. 1 
Integrators S4 PRK64 
r= 11 1/22” 1/53” 
r = 88 1/48” 2/24" 
3.2 ”轨道 动力 学 
图 1 中 的 两 个 轨道 对 同一 算法 的 哈密 顿 误差 的 


长 期 表现 不 同 是 由 于 两 轨道 有 不 同 的 动力 学 性 质 
实际 上 , 初始 半径 r = 11 对 应 有 序 轨道 , 即 轨 道 是 
拟 周 期 的 ; 而 初始 半径 r = 88 对 应 混沌 轨道 , 即 轨 
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运动 


道 演化 对 初始 条 件 微小 改变 有 指数 式 的 敏感 依赖 


性 . 
轨道 的 有 序 和 混沌 动力 学 属性 可 借助 轨道 
的 3 维 空间 构 型 和 平面 投影 图 来 粗略 地 观察 . 选 


取 初 始 轨道 半径 > = 25, LT AE BEE = 0.995 和 
角 动 量 L = 4. 图 2 (a)-(c) 画 出 了 3 个 不 同 磁场 参 
数 8 的 3 维 空间 轨道 构 型 ， 图 中 红色 为 轨道 在 赤道 
Ho = 7/2 上 的 投影 在 图 2 (a) 中 ，B6 = 0 显示 试 
验 粒 子 在 史 瓦 西 黑 洞 时 空中 的 运动 轨道 限于 赤道 
平面 上 ， 并 且 围 绕 黑 洞 做 圆 轨 道 周 期 运动 . 在 图 2 
hh，8 = 1.2x10-4 对 应 3 维 空间 轨道 不 在 一 个 
平面 上 ， 其 投影 是 一 组 相交 环 曲线 . 在 图 2 (c) 中 ， 
B = 0.01 对 应 3 维 空间 轨道 也 不 在 一 个 平面 上 , 其 投 
影 由 一 些 互相 交错 环 曲线 组 成 . 图 2 (c) 与 图 2 (b) 相 
比 , 轨迹 更 杂乱 、 更 没有 规则 . 这 两 个 运动 轨道 是 
混沌 还 是 有 序 难 以 判断 . 轨道 有 序 或 混沌 的 准确 判 
定 还 需要 借助 其 他 方法 
3.2.1 混沌 指标 

判断 带电 粒子 在 史 瓦 西 黑洞 周围 运动 的 轨道 

和 否 混沌 有 许多 方法 . 例如 庞 加 莱 截 面 法 、 李 雅 普 
诺 夫 指数 和 快速 李 雅 普 诺 夫 指标 等 . 

庞 加 莱 截 面 适 合 研究 2 自由 度 的 4 维 保守 系统 ， 
可 用 于 判断 系统 的 运动 状态 . 若 以 0 = x/2 为 曲面 ， 
上 方 为 9 > T/2， 下 方 为 6 < T/2， 可 以 按照 线性 揪 
值 来 求 曲面 上 的 点 . 当 曲面 上 只 得 到 几 个 点 或 者 这 
些 点 构成 一 封闭 曲线 时 , 系统 做 拟 周 期 运动 ; 当 截 
面 上 得 到 的 点 是 杂乱 无 章 随机 分 布 在 一 个 区 域 时 
系统 做 混沌 运动 . 

李 雅 普 诺 夫 指 数 (Lyapunov indicator) 是 衡量 
邻近 轨道 随时 间 平 均 指 数 分 离 比 的 指标 , 能 够 反 
映 轨 道 混 沌 的 强度 . 最 大 李 雅 普 诺 夫 指 数 判断 混沌 
较为 准确 , 计算 方法 主要 有 变 分 法 和 两 粒子 法 21. 
本 文 使 用 变 分 法 计算 : 


地 


T 


IS(7)| 
[S(0)| ° 
别 表示 7 时 刻 与 初始 时 刻 的 切 
个 正 值 , 说 明 有 界 轨 道 处 在 
序 轨道 处 于 有 序 状态 . 


ye (17) 


In 
THO T 


: S(r) 和 S(0) 分 
值 趋 于 稳定 到 
; 入 值 趋 于 0, 说 明 有 


ra 


= 出 
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ii 
状态 


ae 
Si 
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和 是 定义 的 一 个 值 , 用 来 表示 轨道 的 混沌 状态 . 
快速 李 雅 普 诺 夫 指标 (Fast Lyapunov indica- 

tor, FLI) 与 李 雅 普 诺 夫 指 数 相 比 , 可 以 更 快 地 看 到 

轨道 的 混沌 , 得 到 混沌 指数 , 其 定义 为 59l: 


FLI = lg|S(7)|. (18) 


AFLI HAKR, MA FL; A 


FLI 随 时 间 增 长 非常 慢 , 则 有 界 轨道 有 序 . 
(a) B=0 


到 2 ”3 维 轨道 构 型 及 其 在 赤道 面 9 = /2 上 的 投影 


Fig.2 Three-dimensional trajectories and their projections 


at the plane 0 = 7/2 


3.2.2 ”粒子 运动 轨道 的 动力 学 

图 3 为 赤道 面 9 = r/2 上 的 庞 加 莱 截 面 . 图 3 (a) 
为 图 2 (b) 中 粒子 轨道 在 赤道 面 9 = r/2 上 的 庞 加 莱 
AHA, 截面 上 的 轨道 是 闭合 的 环 轨道 ， 说明 图 2 
(b) 中 粒子 做 有 序 的 拟 周期 运动 . 图 3 (b) 为 图 2 (c) 
中 粒子 轨迹 的 庞 加 莱 截 面 图 截面 上 的 轨道 是 随 
机 离散 的 点 分 布 在 一 个 面 区 域 上 ， 表 明 粒 子 做 混 
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沌 运动 . 显然 , 相 比 3 维 空间 构 型 和 平面 投影 图 , 庞 
加 莱 截 面 区 分 轨道 有 序 或 混沌 要 更 清楚 、 准 确 . 带 
BE 试验 粒子 在 史 瓦 西 黑洞 周围 运动 时 , 磁场 破坏 了 
第 4 个 运动 积分 的 存在 , 导致 系统 不 可 积 , 是 粒子 运 
动 混沌 的 根本 原因 . 


P, 


0 50 100 150 200 


10 20 30 40 50 


3 粒子 轨道 在 赤道 面 9 = x/2 上 的 庞 加 莱 截 面 . 图 (a): 图 2 (b) 中 
粒子 的 轨道 . 图 (b): 图 2 (c) 中 粒子 的 轨道 . 


Fig.3 Poincaré cross-section of the particle orbit at the 
plane 0 = 7/2. Panel (a): the particle trajectory is the same 
as that in Fig. 2 (b). Panel (b): the particle trajectory is the 

same as that in Fig. 2 (c). 


图 4 为 进一步 改变 磁场 参数 8 值 时 粒子 轨道 演 
化 情况 . 图 4 (a)-(c) 为 粒子 轨道 运动 的 庞 加 莱 截 面 
图 , 画 了 3 条 轨道 : 轨道 1 (红色 轨道 r = 40), 轨道 2 
( 蓝 色 轨 道 r = 65) 和 轨道 3 (黑色 轨道 r = 88), HR 
轨道 为 绿色 轨道 r = 11， 紫 色 轨 道 r = 20. 固定 粒 
子 能 量 媚 = 0.995 与 角 动 量 L = 4.6 不 变 ， 但 磁场 
参数 在 3 个 子 图 中 依次 为 8 = 1.2x10-4、B = 5.6x 
1074, 6=8.9x10-*. 图 4 (a) 中 6 =1.2x10-4, 所 
有 轨道 都 在 做 有 序 的 拟 周 期 运动 , 因此 截面 上 的 相 
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轨道 都 是 闭合 的 环 轨道 . 图 4 (d) 中 的 最 大 李 雅 普 诺 。 ”表明 此 轨道 为 混沌 轨道 . 在 图 4 (h) 中 轨道 2 快速 李 
夫 指 数 也 同样 揭示 图 4 (a) 中 所 有 轨道 都 是 有 序 的 ， ” 雅 普 诺 夫 指标 随时 间 指 数 增长 ,说明 轨 道 ? 是 混沌 


因 


为 当 积分 时 间 达 


还 不 能 达到 稳定 值 ， 


JJF 
小 于 3， 


DE 


到 107 后 ， 这些 李 雅 普 详 夫 指数 的 . 因此 , 3 种 方法 都 判定 图 4 (b) 中 轨道 1、3 为 有 序 
都 有 趋 于 0 的 趋势 , 表明 轨道 的 轨道 , 轨道 2 为 混沌 轨道 . 图 4 (c) 中 轨道 1 与 轨道 3 为 


生 . 图 4 (g) 中 所 有 轨道 的 积分 到 最 后 FLI 值 都 混沌 轨道 , 轨道 2 形成 5 个 有 序 岛 屿 . 轨道 1、3 的 混 
没有 随时 间 指 数 增长 的 趋势 , 属于 有 序 轨 沌 性 与 轨道 2 的 有 序 性 质 也 被 图 4 (f)、(i) 中 的 最 大 


道 的 特征 . 图 4 (b) 中 6 = 5.6x10~4, 粒子 轨道 2 呈现 李 雅 普 诺 夫 指 数 和 快速 李 雅 普 诺 夫 指 标 所 证 实 . 故 
弱 混 沌 性 , 轨道 1 依 


图 4 (e)、 


旧 做 拟 周期 运动 , 轨道 3 亦 如 此 . 得 出 结论 : 随 着 正 磁场 参数 6 的 增 大 , 亦 即 B62 增 大 ， 


(hp) 中 的 李 雅 普 诺 夫 指 数 和 快速 李 雅 普 诺 。 ”粒子 运动 轨道 从 有 序 走 上 混沌 , 并 且 混 沌 程度 也 随 


夫 指 标 分 别 显示 轨 
图 4 (e) 中 的 李 雅 普 诺 夫 指 数 趋 于 一 个 稳定 的 正 值 ， 力 表现 为 引力 效果 变 大 的 缘故 . 


(a) B=1.2x104 


0 30 60 90 120 150 180 210 


道 1、3 为 有 序 轨道 ,而 轨道 2 在 之 增强 . 这 是 因为 带电 粒子 受到 来 自 磁场 的 洛 伦 兹 


7 


(g) B=1.2x104 


Ig |A| 
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(h) B=5.6x104 
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、 最 大 李 雅 普 诺 夫 指 数 、FLI 随 磁 参 数 的 变化 . 图 (a)-(c): 3 条 研究 轨道 半径 (> = 40,r = 65,r = 88) 的 庞 加 莱 截 面 图 , 粒 


0.995, 角 动 量 L = 4.6; 图 (d)-(f) 为 图 (a)-(c) 的 最 大 李 雅 普 诺 夫 指 数 图 ; 图 (g)-(i) 为 图 (a)-(c) 的 FLI. 


Fig.4 The variation of the Poincaré cross-section, maximum Lyapunov indicator, and fast Lyapunov indicator with magnetic 


parameters. Panels 


a)—(c): Poincaré cross-sections of three studied orbital radii (r = 40,7 = 65, r = 88). Particle energy 


E = 0.995, and angular momentum L = 4.6; panels (d)—(f) are the maximum Lyapunov indicator graphs of panels (a)—(c), 


respectively; panels (g)—(i) are FLI of panels (a)—(c), respectively. 


图 5 中 磁 参 数 8 为 负 值 时 , 粒子 轨道 的 庞 加 羔 截 都 是 拟 周期 有 序 轨道 . 因此 , 随 着 负 磁 场 参数 6 的 绝 


面 随 磁 参 数 的 变化 . 考虑 5 条 轨道 , 初始 半径 依次 对 值 减 小 , 粒子 的 轨道 从 混沌 状态 变 为 有 序 状 态 . 
Ar=1l, r=20. r=40. r=65, r= 88, 颜色 也 就 是 说 , 46 < 0 时 , 8 值 的 增加 , 亦 即 82 变 小 , 会 


is 


图 5 (b 


意义 同 图 4. 图 5 (a) 中 有 两 条 轨道 处 于 强 混 沌 状态 ; 减弱 轨道 的 混沌 程度 . 这 是 由 于 82? 变 小 使 得 洛 伦 兹 
) 中 也 有 两 条 轨道 处 于 混沌 状态 , 但 属于 弱 混 力作 为 引力 效果 也 减 小 . 总 之 , 从 图 4、5 可 以 得 出 


沌 轨道 ; 图 5 (c) 中 只 出 现 一 条 混沌 轨道 , Ab og 无 论 磁 场 参数 6 正 负 , 8 的 绝对 值 变 大 会 加 剧 轨道 的 
混沌 状态 ; 图 5 (d) 中 没有 混沌 轨道 出 现 , 所 有 轨道 混沌 程度 . 


J 


chinaXiv:202308.00109v1 


Chinax iv ERAT 


64 卷 天 文 学 报 4 期 


La -3 
©) B=-1.17x10° | 


Tie, 


“0 20 40 60 80 100 120 0 20 40 60 80 100 120 140 160 
r r 


5 BEBE ATEN, 轨道 的 庞 加 莱 截 面 随 56 的 变化 . 粒子 能 量 妃 = 0.995 与 角 动 量 L = 4.6. 


Fig.5 The variation of the Poincaré cross-section with 8 when the magnetic parameter 8 is negative. Particle energy E = 0.995 


and angular momentum L = 4.6. 


P 


6 为 改变 粒子 能 量 五 的 庞 加 莱 截 面 ， 颜 色 意 致 . 此 外 , 图 6 (b)、(e) 或 图 6 (c)、(f) 再 次 支持 图 4 关 
义 同 图 4. 图 6 (a)-(c) 为 固定 磁 参 数 6 = 8.9x10-4 与 于 正 磁场 参数 变 大 引起 混沌 加 剧 的 结论 . 

角 动 量 L = 4.6， 粒 子 能 量 分 别 为 EB =0.99. E = 图 7 绘制 了 3 个 不 同 粒 子 角 动量 L 的 轨道 演化 行 
0.9925 E = 0.995 的 相 空 间 结 构 情况 .图 6 (a) 展 为 , 颜色 意义 同 图 4. 图 7 (a)-(c) 为 固定 磁场 参数 6 
示 3 条 轨道 都 为 有 序 轨道 . 在 图 6 (b) 中 轨道 3 由 = 8.9x10-* ALF ABLE = 0.995， 粒 子 角 动 量 分 
箭头 型 轨道 与 三 角形 轨道 组 成 ， 其 连接 处 对 别 为 L = 4.21, L= 4.65 L = 4.88 的 庞 加 莱 截 面 情 
应 双 曲 不 动 点 ， 使 得 这 条 轨道 处 于 弱 混 沌 状 形 . 图 7 (a) 显 示 3 条 轨道 (r = 40, r= 65, r = 88) 都 
态 ; 轨道 1 与 轨道 2 为 有 序 轨道 . 图 6 (cj 出 现 2 条 是 混沌 的 , 且 都 为 强 混沌 轨道 . 图 7 (d) 为 图 7 (a) 的 
强 混沌 轨道 ， 即 轨道 1 与 轨道 3， 轨 道 2 由 5 个 岛屿 最 大 李 雅 普 诺 夫 指 数 图 , 图 中 再 次 验证 3 条 轨道 为 
组 成 ， 属 于 有 序 轨道 ， 图 6 (d)-(f) 为 固定 磁 参 混沌 轨道 , 因为 曲线 都 不 趋 于 0, 而 是 趋 于 一 个 稳 
数 B = 1.5x10 习 与 角 动 量 L = 4.6， 粒 子 能 量 分 别 定 正 值 . 图 7 (g) 为 图 7 (a) 的 FLI 图 , 图 中 3 条 轨道 的 
AE =0.99. E = 0.9925 E = 0.995 的 截面 情况 . FLI 值 都 呈 指 数 增长 , 所 以 3 条 轨道 的 混沌 程度 为 强 
在 图 6 (d) 中 轨道 1 与 轨道 2 都 为 马鞍 形状 , 轨道 3 为 混沌 . 图 7 (b) 与 图 4 (c) 的 参数 条 件 相 同 ,轨道 1 与 
箭头 型 轨道 ,3 条 轨道 都 为 有 序 轨道 . 图 6 (e) 中 轨 轨道 3 为 混沌 轨道 , 轨道 2 为 有 序 轨道 . 图 7 (e) 与 图 7 
道 1 与 轨道 2 明显 为 混沌 轨道 , 轨道 3 是 一 个 有 序 环 (h) 分 别 为 图 7 (b) 的 最 大 李 雅 普 诺 夫 指 数 与 快速 李 
轨道 . 图 6 (f) 中 所 有 轨道 都 为 混沌 轨道 . 于 是 ， 从 雅 普 诺 夫 指标 图 . 图 7 (c) 中 轨道 1 为 有 序 轨道 ， 轨 
图 6 中 可 以 得 出 结论 为 粒子 能 量 适当 增加 会 加 剧 轨 道 3 为 混沌 轨道 , 但 轨道 2 为 多 个 岛屿 型 的 有 序 轨道 . 
道 混沌 程度 , 这 是 基于 能 量 增加 会 引起 引力 增强 所 在 图 7 (人 分、 人 的 最 大 李 雅 普 诺 夫 指数 与 快速 李 雅 
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普 诺 夫 指标 图 中 可 以 确定 轨道 2 为 有 序 轨道 . 综 上 BUNA SITE LIE, 点 2 为 弱 混 沌 轨道 ,点 3 表现 为 
所 述 , 得 出 结论 为 : 随 着 粒子 角 动 量 增加 , 轨道 从 弱 有 序 轨 道 . 从 图 8 (a)、(b) 中 得 出 结论 也 能 与 前 面 
混沌 演变 为 有 序 行为 , 从 强 混沌 演变 为 弱 混 沌 . 粒 图 5 的 对 应 , 即 粒子 能 量 B 适 当 增 大 可 以 增强 轨道 
子 角 动量 增 大 减弱 轨道 的 混沌 性 是 源 于 角 动 量 表 混沌 程度 . 图 8 (c)、(d) 绘 制 了 固定 磁场 参数 8 = 
现 为 斥 为 效果 从 而 削弱 引 为 效应. 8.9x10-4 与 能 量 B = 0.995, 而 角 动 量 L 从 3.8 到 5 的 

正如 文献 [22] 所 述 , 除了 用 庞 加 羔 截 面 洞察 动 ”FILI 图 ， 图 8 (c) 中 初始 轨道 半径 x = 11 时 ， 粒 子 
ARERR AES SEER, FLT gin < 41m IRIN, MERTA 
E TARE 8 ee PE EES eg rey OL AR AB MRE AT PE HE 
LETS 2 RUN TORAR, ATR PARET T 点 2 为 有 序 轨道 图 8 (d) 中 初始 轨道 半径 为 r = 88, 


轨道 . 图 8 (a)、(b) 绘 制 了 固定 磁场 参数 8 = 89x 7 a 
0“ 5ABEL = 46, eeeMoosMrLg, WHEAT — RERAN APAS 


e 


= 


= 


每 个 FLI 值 取 自 积分 时 间 到 r = 107 的 结果 ; FLI = 点 2 所 示 的 庞 加 莱 截 面 缘 为 弱 混 沌 轨道 . 从 图 8 (c)、 
3 为 轨道 有 序 与 混沌 的 分 界线 ， 即 FLI > 3 指示 有 (d) 中 得 出 结论 也 能 与 前 面 图 6 的 对 应 , 即 粒 子 角 动 


界 轨 道 的 混沌 性 ,而 PLI 不 超过 3 表明 有 界 轨道 的 。 量 了 变 入 可 以 减 小 轨 关 混沌 程度 . 图 8 (e). (Pat 
JFE. 图 8 (a) 中 初始 轨道 半径 x = 11, EDF 了 固定 粒子 角 动 量 L = 4.6 与 能 量 B = 0.995, 磁场 
0.9947 为 有 序 轨道 , 巨大 于 0.9957 时 混沌 开始 ; 图 中 。 ”参数 8 从 0 跨度 到 1.6 x 10 的 FLI 图 . SCR TERY Ih 
点 3 所 示 的 庞 加 莱 截 面 为 弱 混 沌 轨道 ,点 1 与 点 2 为 界 6 值 依次 为 9.01 x 10-4、6.92 x 10-4. 从 图 8 (e), 
I 序 轨道 . 图 8 (b) 中 初始 轨道 半径 > = 88, IKE = (f) 中 得 出 结论 与 前 面 图 3、4 的 对 应 , 即 随 着 正 磁场 
0.9928 大 的 能 量 诱发 混沌 ; 图 中 点 1 所 示 的 庞 加 莱 参数 6 的 增加 会 增强 轨道 混沌 程度 . 
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A6 改变 粒子 能 量 刀 的 庞 加 莱 截 面 图 , 粒子 角 动 量 L = 4.6. 


Fig.6 The variation of the Poincaré cross-section with the particle energy E, particle angular momentum L = 4.6. 
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7 ”轨道 的 庞 加 莱 截 面 、 最 大 李 雅 普 诺 夫 指 数 、FLI 随 试验 粒子 角 动量 工 的 变化 . WEBB = 8.9x1074, 粒子 能 量 B = 0.995. 
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Fig.7 The variation of the Poincaré cross-section, maximum Lyapunov indicator and fast Lyapunov indicator with the angular 


momentum of experimental particle L. Magnetic parameter 8 = 8.9x 1074 and particle energy E = 0.995. 
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图 8 快速 李 雅 普 诺 夫 指 标 FLI 与 参数 的 依赖 关系 . 每 个 FLI 值 取 自 积分 时 间 r = 107 后 的 结果 , 有 序 与 混沌 的 FLII 临 界 值 为 3. 图 (a)、 
(b) HL = 4.6, 8 = 8.9x10-4. A(c), (d)#E = 0.995, 8 = 8.9x107*. 图 (e)、(f) 中 LL = 4.6. E = 0.995. 


Fig.8 The dependence between fast Lyapunov indicator FLIs and parameters. Each of the FLIs are obtained after the 
integration 7 = 10", and the FLI value from order to chaos is 3. In panels (a) and (b), L = 4.6 and 6 = 8.9x 10741; in panels (c) 
and (d), E = 0.995 and 6 = 8.9x 107f; in panels (e) and (f), L = 4.6 and E = 0.995. 
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间 对 应 的 FLI 分 布 . 图 9 (a) 是 2 
应 的 FLI 分 布 . 初始 轨道 半径 


s 


Veli Ae EE MA E 


0.995772 AF aR AEE, 并 
的 增 大 而 增 大 , 增 大 到 一 定 值 


后 混沌 开始 减弱 ， 呈 现在 图 中 颜色 由 蓝 色 依次 过 


渡 到 红色 , 再 由 红色 


EE 


的 一 维 图 是 一 致 的 . 


可 以 减弱 轨 i 


BHE = 0.995 时 显示 


道 半径 > = 88 时 扫描 参数 已 和 厂 的 2 维 分 布 图 情形 ， 


外 0.996 


= 0.9957 到 已 = 


转换 到 蓝 色 ， 而 混沌 区 域 大 至 
1 之 间 ， 这 个 结果 与 图 8 (a) 
还 可 以 看 出 增 大 粒子 角 动量 


道 的 混沌 程度 , 这 与 图 8 (c) 的 一 维 图 ， 


的 情况 相同 . 图 9 (c) 为 初始 轨 
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更 清晰 地 展示 了 粒子 角 动 量 增 大 减弱 了 轨道 的 混 
沌 程度 ， 这 与 图 8 (b)、(d) 显 示 的 结果 相同 . 扫描 
图 9 (b) 中 参数 已 和 8 的 2 维 分 布 图 显示 在 初始 轨道 
半径 r = 11 时 , 能 量变 大 , 故 粒 子 受到 的 引力 变 大 
会 加 剧 混沌 程度 ; 但 能 量 大 到 接近 1 时 , 轨道 有 不 稳 
定 趋势 , 相应 减弱 混沌 程度 ; 而 当 6 值 增 大 时 , 轨道 
由 有 序 变 为 混沌 , 且 FLI 值 逐渐 变 大 , 即 混沌 程度 逐 
渐 加 强 . 图 9 (d) 为 初始 轨道 半径 > = 88 时 扫描 参 
Soe Al, 可 以 得 到 与 图 9 (b) 相 同 的 结 

)、(d) 与 图 8 (e)、(f) 扫 描 的 一 维 图 显示 
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B(10"%) 
图 9 ”2 维 参数 空间 对 应 的 FLI 分 布 . 积分 时 间 均 为 r = 107; 图 (a)、(c): 8 = 8.9x1074; 图 (b)、(d): L = 4.6. 


Fig.9 The corresponding FLI distribution for two-dimensional parameter spaces. The integral time is 7 = 107 for all plots; 
panels (a), (c): 8 = 8.9x 107‘; panels (b), (d): L = 4.6. 
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总 结 与 展望 
辛 算法 因 具 有 保 结构 性 质 , 可 为 天 体 长 期 定性 
化 研究 提供 可 靠 的 数值 结果 . 在 计算 效率 方面 


MEL 
DR 


分 ， 


S EAA I: 
部 分 分 


5J 
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系 
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密 顿 函数 可 以 分 解 


显 式 辛 算法 比 隐 式 辛 算法 好 , 磁化 史 瓦 西 时 空 的 哈 


为 具有 显 式 可 积 解 的 93、4、5 部 
可 以 进行 显 式 辛 算法 的 构建 与 应 用 . 优化 
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的 磁场 下 史 瓦 西 黑洞 周围 带电 粒子 轨道 动力 学 . 

工作 主要 探讨 磁场 参数 5、 粒子 能 量 B 和 和 角 动 


量 L 对 粒子 轨道 动力 学 行为 的 影响 . 磁场 是 哈密 顿 


统 不 可 积 性 
值 的 变 大 可 以 使 得 粒子 运 


和 产生 混沌 的 关键 因素 . 磁场 参数 绝 
动 轨道 由 有 序 变 为 混 


E, 增强 轨道 混沌 程度 . 粒子 能 量 增加 亦 可 加 强 混 


沌 程度 , 但 粒子 角 动 量 增 大 反 


而 对 混沌 有 抑制 作用 
于 许多 复杂 的 相对 


本 工作 所 建立 的 算法 适 


论 时 空 的 哈密 顿 函数 或 时 间 变换 哈密 顿 函数 , 可 以 


FA 


来 探索 黑洞 轨道 动力 学 性 质 , 还 可 以 用 于 黑洞 阴 


影 数值 模拟 研究 . 
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Chaotic Motion of Charged Test Particles in a Magnetized 
Schwarzschild Black Hole 


ZHOU Na-ying' ZHANG Hong-xing' SUN Xin? LIU Wen-fangl LI Dan? 


(1 Center of Application and Research of Computational Physics, School of Mathematics, Physics and Statistics, Shanghai 
University of Engineering Science, Shanghai 201620) 
(2 Guangxi Key Laboratory for Relativistic Astrophysics, Guangxi University, Nanning 530004) 


Asstrract The Hamiltonian describing the motion of charged particles around the Schwarzschild black 
hole immersed in an external magnetic field is nonintegrable. Such relativistic Hamiltonian systems do not 
have two splitting parts with analytical solutions as explicit functions of time. This leads to the difficulty 
in the construction and application of explicit symplectic algorithms to the relativistic systems. Recently, 
Chinese scholars have published a series of works in the Astrophysical Journal, where explicit symplectic 
methods are successfully designed for these relativistic Hamiltonians split into three or more explicit 
integrable parts. There are two questions of whether the numbers of splitting these Hamiltonians affect the 
numerical accuracy and which of the explicit symplectic integrators shows the best performance. Our latest 
work in the Astrophysical Journal answered the two questions, and shows that the fourth-order optimal 
Partitioned-Runge-Kutta (PRK¢4) explicit symplectic algorithms with the three-part splitting method as 
the least number of splitting these Hamiltonians performs the best accuracy. This paper applies such an 
integrator to obtain Poincaré cross-section, maximum Lyapunov indicators and fast Lyapunov indicators 
(FLIs), which distinguish between the regular and chaotic dynamical properties of charged particles moving 
near the magnetized Schwarzschild black hole. For given specific values of the particle energy and angular 
momentum, a small magnetic field does not induce chaos, whereas a large positive magnetic field parameter 
easily causes the occurrence of chaos. The strength of chaos increases with the magnetic field increasing. 
Chaos is also strengthened as the particle energy increases. However, it is weakened when a negative 
magnetic field parameter and the particle angular momentum increase. 


Key words celestial mechanics, black hole, magnetic field, chaotic, integrators 
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